// atmos_profile.cpp

#include <cmath>
#include <algorithm>

#include <quigon/basic/atmos_profile.h>
#include <quigon.h>

using namespace std;

namespace quigon
{
Atmos_Profile load_atmos_profile(
        const string& filename,
        size_t size,
        double bin_size,
        double start_value
)
{
	try
	{
		Serie<double, double> raw_serie(filename);
		return serie_to_profile<Atmos_Profile> (
		               raw_serie,
		               size,
		               bin_size,
		               start_value
		       );
	}
	catch (...)
	{
		func_err("load_atmos_profile", filename);
		throw;
	}
}

Atmos_Profile cal_beta_m(
        const Atmos_Profile& t_profile,
        const Atmos_Profile& p_profile,
        double lambda,
        double rho_n
)
{
	try
	{
		if (
		        ! float_equal(t_profile["bin_size"], p_profile["bin_size"]) ||
		        ! float_equal(t_profile["start_value"], p_profile["start_value"])
		)
			throw invalid_argument("error: different attr profiles used to make beta_m");

		size_t final_size = min(t_profile.size(), p_profile.size()) ;
		Atmos_Profile result(final_size);

#ifdef VERBOSE
		cout << "  Calculating BetaM..." << endl;
#endif
		//double n0 = 1.000272599; //(n0-1)*10e8=6432.8+2949810/(146-lambda^-2)+25540/(41-lambda^-2);
		//Edlen 1953, I doubt about it though
		double n0_minus_1 =
		        (
		                6432.8 + 2949810.0 / (146.0 - pow(lambda, -2.0))
		                + 25540.0 / (41.0 - pow(lambda, -2.0))
		        )
		        * 1.0e-8;
		//notice: in this equation, lambda is used as nm.

		//test
#ifdef VERBOSE
		cout << "n0 - 1 = " << n0_minus_1 << endl;
#endif
		//BetaM(z)=PI^2/(N0^2*lambda^4) * (n0^2-1)^2 * (6+3*rho_n)/(6-7*rho_n) * N0 * T0/P0 * P(z)/T(z)
		//and n0^2-1 can be written into 2(n0-1)
		//where N0=2.547e25 , T0=288, P0=1013.25, rho_n=0.0284, lambda=523?
		//and PI^2/N0 * T0/P0 =1.1014e-25
		for (size_t i = 0; i < final_size; ++i)
		{
			//Well, take care of the Units here.
			result[i] =
			        1.1014e-25 * pow((2.0 * n0_minus_1), 2)
			        / pow(lambda * 1.0e-9, 4)
			        * (6.0 + 3.0 * rho_n) / (6.0 - 7.0 * rho_n)
			        * p_profile[i] / t_profile[i] * 1000.0; // notice beta_m is in km^-1 term.

		}
		//and here lambda is used as meter, so there's a 1.0e-9, and result is in km^-1, so multiplies 1000.0.
		return result;
	}
	catch (...)
	{
		func_err("cal_beta_m");
		throw;
	}
}

}


